* fig1-estimate-garch
* 4/17/24
*
* ---------------------------------------------------------------------------
clear
set maxvar 8000
cd "~/Dropbox/Benton-Jordan-Philips/final-JOP-replication"
use "data/simulated_wide.dta", clear

tsset time

gen estly = .
gen estlx = .
gen estx = .
gen estalph = .
gen esthetx = .
gen esthetz = .
gen estomega = .
gen estarch = .
gen estgarch = .
gen estlyse = .
gen estlxse = .
gen estxse = .
gen estalphse = .
gen esthetxse = .
gen esthetzse = .
gen estomegase = .
gen estarchse = .
gen estgarchse = .

* test 1 garch
arch y_comb_2 l.y_comb_2 x_mean l.x_mean, arch(1) garch(1) het(x_het_2 z_het_2) 
mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B

forvalues i = 1(1)480 {
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	capture replace estly = [y_comb_`i']_b[L.] if combination == `i'
	capture replace estlx = [y_comb_`i']_b[x_mean] if combination == `i'
	capture replace estx = [y_comb_`i']_b[L.x_mean] if combination == `i'
	capture replace estalph = [y_comb_`i']_b[_cons] if combination == `i'
	capture replace esthetx = [HET]_b[x_het_`i'] if combination == `i'
	capture replace esthetz = [HET]_b[z_het_`i'] if combination == `i'
	capture replace estomega = [HET]_b[_cons] if combination == `i'
	capture replace estarch = [ARCH]_b[L.] if combination == `i'
	capture replace estgarch = [ARCH]_b[L.garch] if combination == `i'
	capture replace estlyse = [y_comb_`i']_se[L.] if combination == `i'
	capture replace estlxse = [y_comb_`i']_se[x_mean] if combination == `i'
	capture replace estxse = [y_comb_`i']_se[L.x_mean] if combination == `i'
	capture replace estalphse = [y_comb_`i']_se[_cons] if combination == `i'
	capture replace esthetxse = [HET]_se[x_het_`i'] if combination == `i'
	capture replace esthetzse = [HET]_se[z_het_`i'] if combination == `i'
	capture replace estomegase = [HET]_se[_cons] if combination == `i'
	capture replace estarchse = [ARCH]_se[L.] if combination == `i'
	capture replace estgarchse = [ARCH]_se[L.garch] if combination == `i'
	display `i'
}

save "data/garch-estimates-Stata.dta", replace


* make directories to store datasets:
cap mkdir "data/figure1-data"
	cap mkdir "data/figure1-data/hetxatzero"
		cap mkdir "data/figure1-data/hetxatzero/hetzatzero"
			cap mkdir "data/figure1-data/hetxatzero/hetzatzero/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatzero/shockpos1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatzero/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatzero/shockneg1"
		cap mkdir "data/figure1-data/hetxatzero/hetzatpos1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatpos1/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatpos1/shockpos1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatpos1/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatpos1/shockneg1"
		cap mkdir "data/figure1-data/hetxatzero/hetzatneg1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatneg1/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatneg1/shockpos1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatneg1/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatzero/hetzatneg1/shockneg1"
	
	cap mkdir "data/figure1-data/hetxatpos1"
		cap mkdir "data/figure1-data/hetxatpos1/hetzatzero"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatzero/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatzero/shockpos1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatzero/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatzero/shockneg1"
		cap mkdir "data/figure1-data/hetxatpos1/hetzatpos1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatpos1/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatpos1/shockpos1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatpos1/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatpos1/shockneg1"
		cap mkdir "data/figure1-data/hetxatpos1/hetzatneg1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatneg1/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatneg1/shockpos1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatneg1/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatpos1/hetzatneg1/shockneg1"
	
	cap mkdir "data/figure1-data/hetxatneg1"
		cap mkdir "data/figure1-data/hetxatneg1/hetzatzero"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatzero/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatzero/shockpos1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatzero/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatzero/shockneg1"
		cap mkdir "data/figure1-data/hetxatneg1/hetzatpos1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatpos1/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatpos1/shockpos1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatpos1/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatpos1/shockneg1"
		cap mkdir "data/figure1-data/hetxatneg1/hetzatneg1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatneg1/shockpospoint1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatneg1/shockpos1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatneg1/shocknegpoint1"
			cap mkdir "data/figure1-data/hetxatneg1/hetzatneg1/shockneg1"














* 36 combinations
*  1-4 . hetx at -1; hetz at -1; shock [1]  -1, [2]  -0.1, [3]  0.1, [4]  1
*  5-8 . hetx at -1; hetz at  0; shock [5]  -1, [6]  -0.1, [7]  0.1, [8]  1
*  9-12. hetx at -1; hetz at +1; shock [9]  -1, [10] -0.1, [11] 0.1, [12] 1

*  1-4 . hetx at  0; hetz at -1; shock [13] -1, [14] -0.1, [15] 0.1, [16] 1
*  5-8 . hetx at  0; hetz at  0; shock [17] -1, [18] -0.1, [19] 0.1, [20] 1
*  9-12. hetx at  0; hetz at +1; shock [21] -1, [22] -0.1, [23] 0.1, [24] 1

*  1-4 . hetx at +1; hetz at -1; shock [25] -1, [26] -0.1, [27] 0.1, [28] 1
*  5-8 . hetx at +1; hetz at  0; shock [29] -1, [30] -0.1, [31] 0.1, [32] 1
*  9-12. hetx at +1; hetz at +1; shock [33] -1, [34] -0.1, [35] 0.1, [36] 1







set seed 50

* Do the shocking!. hetx -1; hetz -1; shock -1 (Combination #1)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = z_het_i
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + -1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatneg1/shockneg1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz -1; shock -0.1 (Combination #2)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + -0.1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatneg1/shocknegpoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz -1; shock 0.1 (Combination #3)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + 0.1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatneg1/shockpospoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz -1; shock 1 (Combination #4)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + 1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatneg1/shockpos1/comb_`i'.dta", replace
	}
}



set seed 50

* Do the shocking!. hetx -1; hetz 0; shock -1 (Combination #5)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + -1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatzero/shockneg1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz 0; shock -0.1 (Combination #6)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + -0.1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatzero/shocknegpoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz 0; shock 0.1 (Combination #7)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + 0.1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatzero/shockpospoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz 0; shock 1 (Combination #8)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + 1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatzero/shockpos1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz 1; shock -1 (Combination #9)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + -1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatpos1/shockneg1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz 1; shock -0.1 (Combination #10)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + -0.1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatpos1/shocknegpoint1/comb_`i'.dta", replace
	}
}



set seed 50

* Do the shocking!. hetx -1; hetz 1; shock 0.1 (Combination #11)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + 0.1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatpos1/shockpospoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx -1; hetz 1; shock 1 (Combination #12)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(-1 + 1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(-1) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatneg1/hetzatpos1/shockpos1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz -1; shock -1 (Combination #13)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + -1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatneg1/shockneg1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz -1; shock -0.1 (Combination #14)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + -0.1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatneg1/shocknegpoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz -1; shock 0.1 (Combination #15)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + 0.1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatneg1/shockpospoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz -1; shock 1 (Combination #16)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + 1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatneg1/shockpos1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz 0; shock -1 (Combination #17)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + -1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatzero/shockneg1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz 0; shock -0.1 (Combination #18)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + -0.1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatzero/shocknegpoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz 0; shock 0.1 (Combination #19)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + 0.1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatzero/shockpospoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz 0; shock 1 (Combination #20)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + 1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatzero/shockpos1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz 1; shock -1 (Combination #21)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + -1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatpos1/shockneg1/comb_`i'.dta", replace
	}
}



set seed 50

* Do the shocking!. hetx 0; hetz 1; shock -0.1 (Combination #22)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + -0.1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatpos1/shocknegpoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz 1; shock 0.1 (Combination #23)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + 0.1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatpos1/shockpospoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 0; hetz 1; shock 1 (Combination #24)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(0 + 1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(0) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatzero/hetzatpos1/shockpos1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz -1; shock -1 (Combination #25)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + -1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatneg1/shockneg1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz -1; shock -0.1 (Combination #26)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + -0.1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatneg1/shocknegpoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz -1; shock 0.1 (Combination #27)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + 0.1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatneg1/shockpospoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz -1; shock 1 (Combination #28)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + 1) + b6*(-1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(-1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatneg1/shockpos1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz 0; shock -1 (Combination #29)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + -1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatzero/shockneg1/comb_`i'.dta", replace
	}
}



set seed 50

* Do the shocking!. hetx 1; hetz 0; shock -0.1 (Combination #30)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + -0.1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatzero/shocknegpoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz 0; shock 0.1 (Combination #31)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + 0.1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatzero/shockpospoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz 0; shock 1 (Combination #32)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + 1) + b6*(0) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(0) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatzero/shockpos1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz 1; shock -1 (Combination #33)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + -1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatpos1/shockneg1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz 1; shock -0.1 (Combination #34)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + -0.1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatpos1/shocknegpoint1/comb_`i'.dta", replace
	}
}


set seed 50

* Do the shocking!. hetx 1; hetz 1; shock 0.1 (Combination #35)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + 0.1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatpos1/shockpospoint1/comb_`i'.dta", replace
	}
}



set seed 50

* Do the shocking!. hetx 1; hetz 1; shock 1 (Combination #36)
forvalues i = 1(1)480 {
	clear all
	set maxvar 8000
	use "data/simulated_wide.dta", clear
	tsset time 
	capture qui arch y_comb_`i' l.y_comb_`i' x_mean l.x_mean, arch(1) garch(1) het(x_het_`i' z_het_`i') nolog
	if !_rc {
		* su x_het_1, meanonly
		* loc x_het_themean = r(mean)    				
		* I can't get this to work, so I'll manually set it to 0 (which is the ``truth'' mean and the average value of the averages
		mat B = e(b)
		mat VCV = e(V)
		loc names: colnames VCV
		di "`names'"
		mat list B
		clear	
		corr2data b1 b2 b3 b4 b5 b6 b7 b8 b9, n(1000) means(B) cov(VCV)	
		* 1 = l.y_comb_i
		* 2 = x_mean
		* 3 = l.x_mean
		* 4 = _cons (mean)
		* 5 = x_het_i
		* 6 = het_z
		* 7 = _cons (het)
		* 8 = arch
		* 9 = garch
	
		* First and t to 30, keep at mean (b5 het_x = 0) **************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation
		gen sigma2_0 = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*0
	
		forv j = 1/30 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*r(mean) 
		}

		* Shock at t=31, shock x_het 
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation and b5 to shock size ****************
		su sigma2_30, meanonly
		gen sigma2_31 = exp(b7 + b5*(1 + 1) + b6*(1) + b8*0) + b9*r(mean)

		* Post-shock at t = 32+, keep at mean (b5 het_x = 0)
		**************** CHANGE VALUE ON b5 and b6 to het_x and het_z for simulation ****************
		forv j = 32/50 {
			loc pastj = `j' - 1
			su sigma2_`pastj', meanonly
			gen sigma2_`j' = exp(b7 + b5*(1) + b6*(1) + b8*0) + b9*r(mean)
		}

		* take means and draw CIs:
		forv j = 1/50 {
			_pctile sigma2_`j', p(2.5, 5, 12.5, 87.5, 95, 97.5) 
			gen sigma2_ll95_`j' = r(r1)
			gen sigma2_ll90_`j' = r(r2)
			gen sigma2_ll75_`j' = r(r3)
			gen sigma2_ul75_`j' = r(r4)
			gen sigma2_ul90_`j' = r(r5)
			gen sigma2_ul95_`j' = r(r6)
		}
		* Collapse, reshape, drop 'burnin' values, and plot:
		qui collapse sigma2*
		qui gen count = _n
		qui reshape long sigma2_ sigma2_ll95_ sigma2_ll90_ sigma2_ll75_ sigma2_ul75_ sigma2_ul90_ sigma2_ul95_, j(time) i(count)
		qui keep in 29/50
		qui drop time
		qui gen time = _n - 1
		**************** CHANGE FILE PATH TO CORRECT FOLDER ****************
		qui save "data/figure1-data/hetxatpos1/hetzatpos1/shockpos1/comb_`i'.dta", replace
	}
}




